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Abstract 



In a previous paper, we have described a method to perform Fixed-Node 
Quantum Monte Carlo calculations for lattice fermions. In this paper, we 
present an extension of this method, by which it is possible to find information 
on the properties of the exact ground-state wave function. We give some 
further illustrations of the FNMC and the extended methods. 

PACS numbers: 71.20.Ad, 75.10.Lp, 71.10.-hx 
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I. INTRODUCTION 



Recently, we described a method to perform Quantum Monte Carlo (QMC) calculations 
for lattice fermionJlli. We introduced an effective Hamiltonian, by which the sign problem 
that is usually encountered in such systems, can be avoided. We proved that it leads to an 
upper bound for the ground-state energy. This Fixed-Node Monte Carlo (FNMC) method 
for lattice systems is very much inspired by the similar method, developed by Ceperley et 
a/.ifl, for continuous systems. 

In this paper we want to show that, in analogy to the method for continuum problems, it 
is possible to extend the FNMC method for lattice systems in order to obtain approximate 
information about the true ground state of the system. This nodal relaxation method, as it is 
called by Ceperley et a/.i, can be straighforwardly implemented also on lattice systems. We 
will give a formalism for this method, using the (theoretical) results of the FNMC method 
described in0i. Also, we will give a few more examples of applications of these methods, 
in comparing our results to those presented by other authors, for a few systems frequently 
used in the literature. 

In the next section we will recall the basic formulas of the Green Function Monte Carlo 
(GFMC) formalism, and give our notation and interpretation of the concept of importance 
sampling. After a few summarizing remarks on the Fixed- Node Monte Carlo method, we will 
show in Section |ITl| how the results of this method can be used to calculate properties of the 
true ground state of the system. In Section |rV| we present a few examples of applications 
of both methods, and we conclude this paper by a discussion on the applicability of the 
methods and the need for further investigation of possible wave functions. 

A much more elaborate explanation and discussion of the FNMC method and the ex- 
tended method, information on the implementation of both methods, and some more exam- 
ples of applications, have been given in Refs. ^ and ^ 



II. BASIC FORMULAS OF GREEN FUNCTION MONTE CARLO 

The GFMC method, as we use it, is based on the projection operator 

^= 1 -t(7i:-u;), (1) 

by which the ground state of the Hamiltonian Ti is filtered from a trial state |\&t), which is 
known or can be calculated for each configuration R in the configuration space {R}. The 
parameters r and w are to be chosen appropriately. After applying this operator n times, 
an approximation 

|^("))=jr"|^^) (2) 

of the ground state is obtained, which closely resembles the ground state for large n. The 
energy of this state can be formally calculated as 

E7^(^T|^|i^n) [m=ARmRi-l)] (R^I^t) . . 

E^^{^T\Rn)[m=l{R^\m^-l)]{Ro\'^T) ' 
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where = {Rq, Ri, R2, Rn} is a path in the configuration space. As we will see in the 
actual implementation of this calculation, the wave function has to be known in the end 
point of each path, and therefore {^t\ must be used in this expression in stead of 
This is called a mixed estimate of the energy. 

The expression can usually not be calculated exactly. In the Monte Carlo scheme, 
only a limited number of paths TZ is used for the calculation of this expression. We can 
rewrite (^) as 

(„) ^ E7^igT(i^n)m(7^)p(7^) 

E7^"^(7^)p(7^) ' ^ ' 

where 

^^^^^ - {^^\R) 

is the local energy in the configuration R. The quantity p{TZ) can be interpreted as the 
probability of choosing a complete path TZ, and each path carries a weight or multiplicity 
m{TZ). If one makes sure that p and m satisfy 



p{n)m{n) = {^^\Rn) 

the energy can be calculated as 



\{{R,\F\R,_^) 



u=i 



(^oI^t), (7) 



p(n) _ Y!nET{Rn)m{n) 

where the summation is now restricted to some number of paths, chosen according to the 
probability piTZ)- In the limit of using infinitely many paths, this calculation becomes exact; 
in case of a finite number of paths, a statistical error bar can be defined. 

In the definition of p and m one can include ways to make the sampling of the wave 
function more efficient. This is called importance sampling, and usually refers to the idea 
that one should try to sample more often in regions where the trial wave function has a 
large absolute value. We give a description which includes the possibility of performing 
importance sampling by means of a guiding function \1/g- This function must be positive 
everywhere; often, it is taken to be the absolute value of the trial wave function. We denote 
the probability for choosing a path as 

n 

Pguiding(7^) = PoiRo)Y[p{Ri Ri-l), (9) 
i=l 

Here, the probability for the first configuration of the path to be R is 

Po{R) = {^t\R){R\^t), (10) 
and the probability of each subsequent step of the path is defined by 
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p{E! ^ R) 



{mG\R'){R'\F\R) 



{^!G\R)m{R) 

When using this probabiUty for the paths, one has to take m as 



(11) 



'^guiding ('^) 

where the weight factors m(i?j_i) 

m{R) = y: 

R' 



(^T|i?o) 



.i=l 



(^T|i?n) 
(^G|i?n)' 



can be calculated for each step Ri 
{^g\R'){R'\F\R) {^g\F\R) 



{^g\R) 



(*g|^) 



(12) 

Ri-i separately as 
(13) 



These relations basically define the Monte Carlo procedure. Additionally, one can make use 
of branching^, to avoid that the variance of the calculated result increases too rapidly. 

In the Monte Carlo scheme, one can thus choose the starting configurations of the paths 
according to po, and each subsequent configuration according to the stochastic matrix defined 
by p{R' R). Interpreting the paths as random walks in the configuration space, one 
usually calls the configurations chosen walkers. In our interpretation, the starting set of 
walkers that one uses in the Monte Carlo calculations is distributed according to the square 
of the trial wave function. One can calculate how the walkers are redistributed after n steps, 
considering also the weight they carry: 



piR) = rn{n)pin). 

{n\R^=R} 

This expression can be easily evaluated using the relation 



(14) 



PiR) = {^t\R) E 

{Rf3,Rl,-,Rn-l} 

= (*T|i?)(i?|J^"|^T) 



Yl{R,\F\R,.i) 



i=l 



for p and m: 

(i?o|^T) 



(15) 

(16) 
(17) 



Thus, after a sufficiently large number of steps the walkers are distributed according to the 
product of the trial wave function and the ground-state wave function of the Hamiltonian 
considered. 



III. EXTENSION OF THE FNMC FORMALISM 

The Fixed-Node Monte Carlo method uses in fact the same formalism as the general 
Green Function Monte Carlo, presented in the previous section, but with a modified Hamil- 
tonian TicfT to avoid steos that could introduce a change of sign in the weights or the probabil- 
ities (the sign problemB). The GFMC procedure is straightforwardly applied to this effective 
Hamiltonian, for which a prescription has been given in Refs. |1] and ^ Thus, by the FNMC 
method, information is obtained not on the ground state of the Hamiltonian Ti one is inter- 
ested in, but on the ground state of Heg. It has been proven that the ground-state energy 
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of Tied provides an upper bound for the ground-state energy of H. Of course, it would be 
better to obtain information directly on 7i. As we will see, under certain circumstances this 
is possible, using the FNMC result as input for a new calculation, called nodal relaxation. It 
strongly resembles the method described by Ceperley and Aldem for a continuum problem. 
The method suffers from the sign problem, but it may yield convergent results before the 
fluctuations destroy the accuracy of the calculations. 

Let us rewrite the expression (|T^ to give the distribution of the configurations in the 
FNMC procedure after n steps: 

pS(i?) = (^T|i?)(i?|^S^), (18) 

where {"^l^) is supposed to be a good approximation of the ground state |\I'efr) of Hcs- The 
idea is now that |\E'efr) is sufficiently close to the ground state |\E'o) of H, such that in GFMC 
procedure used on Ti convergence can be obtained within a relatively small number of steps. 
Thus, we modify (H) in the following way: 

|4"))=,F"|vl/eff). (19) 

Then, we also have to modify the expression (^) for the mixed estimate of the energy: 

^(n) ^ (^t|^|4"^) .201 
D \ Inn / D I Z?! D \1 / D liTr _A 

(21) 

Note that we still have to put (^'tI in this expression, as that is the only available information 
we have in each configuration. However, in the starting configurations of the paths, the trial 
wave function has been replaced by the effective ground-state wave function, of which we 
have information through the distribution of the walkers in the FNMC procedure. We 
have to pay further attention to the fact that sign changes can now occur, which have to be 
embedded in the transition probabilities and multiplicities. We can use the following guiding 
procedure, adapted from the regular GFMC scheme described in the previous section: the 
probability to start in a configuration R is given by the fixed-node resulting distribution: 

Po{R)=p'pi{R). (22) 

The transition probabilities are given by 

i^HclR') {R'\F\R) 
piR' ^R)= L oL o ^ (23) 





H\Rn) [Y{U{R^\F\R^-l)] (^C 


l^cff) 




r 


Rn) [Y{U{Ri\F\R.-i)] {Ro\ 





{^>G\R)m{R) 



with the modified weight factor 



_ (^gI^') {R'\f\R) 

m{R)=T. /l J • (24) 



R' 



{^g\R) 
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The absolute value has to be taken to make sure that the transition probabilities remain 
positive, such that the procedure can still be interpreted as a stochastic walk. With these 
expressions replacing (|1^), (0), and ([I3|), the equations (||) and (|1^) can again be used for 
the total probability and the total "weight" of a path, respectively. Note that the weights 
can now be negative, due to the possibly different signs of the wave function in the starting 
and end configurations of the path. 

By this formalism, one can thus try to obtain information directly on the true ground 
state of the Hamiltonian, using the fixed- node result as a starting point. In the following 
section, we will show a few examples of applications of this method. 

IV. APPLICATIONS 

First, we visualize how the FNMC and the nodal-relaxation method work for the Hubbard 
model on a small system, the 2x2x2 cube, with U = 2.5. We have performed exact 
calculations for this syste m0, in order to be able to compare the Monte Carlo results and 
check whether the program has been correctly implemented. In Figures |I| and 0, we show 
calculations for this system at half filling and with zero total spin in the z-direction. In 
Figure |I], the first 250 steps in the Fixed-Node Monte Carlo calculation, using a homogeneous 
mean-field wave function as trial wave function, are shown. As can be seen, after some 100 
steps the energy measured starts fluctuating around the exact ground-state energy of the 
effective Hamiltonian (indicated by the drawn line). The square at the right indicates the 
energy that one obtaines after averaging over several thousands of steps (not shown in 
this figure), with a statistical error which is smaller than the symbol size. In Figure |^, 
two separate calculations of nodal relaxation are shown; one using the same trial wave 
function as in Figure |I| (indicated by H), the other using a mean- field wave function with 
antiferromagnetic ordering imposed on the spins (indicated by AF). The horizontal lines 
indicate the exact ground-state energies of the effective (FN) Hamiltonians using the H and 
the AF trial functions, respectively, and of the true Hamiltonian. The resulting energy after 
each step is indicated by an error bar (the energy measured being the middle value of the 
error bar). As can be seen, the size of the error bars in both calculations first decreases 
with increasing number of steps, and then starts increasing, indicating the existence of the 
sign problem. The calculation marked H seems to converge to the correct value before the 
sign problem arises; the AF calculation does not has not converged at that point. Clearly, 
this is due to the fact that the fixed-node energy found when using the AF wave function 
is significantly higher than the resulting energy in the other fixed-node calculation. This 
nicely demonstrates that the exact result can be reached, but only if the fixed-node result 
is close enough to the true ground state. 

We have performed some more calculations for the Hubbard model. In Table |I|, we 
compare results of calculations of the ground-state energy for a 4 x 4 square with periodic 
boundary conditions (PBC), with 5 spins up and 5 down, and [7 = 4. For this system, exact 
and approximate results are available in the literature. As can be seen, in the case of U = 4, 
the fixed-node method yields an upper bound for the ground-state energy which is quite 
close to the exact value, and nodal relaxation leads to the exact result with an error bar 
of only 0.05% . In Table |I|, we show the results of various calculations of the ground-state 
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energy in a 10 x 10 system (PBC), at half filling and with zero total spin in the z-direction, 
also for [7 = 4. There are no exact results available for this system, but several results of 
mean-field and quantum Monte Carlo calculations can be found in the literature. As can 
be seen, our results are in very good agreement with previous QMC results by Hirsch0 and 
White et aM 



V. DISCUSSION AND OUTLOOK 

First we would like to make a few additional remarks on the calculations we have per- 
formed. As one can see from Table |I|, adding a Gutzwiller factor to the bare mean-field 
wave function greatly improves the (variational) energy obtained from the mean-field cal- 
culation. In the fixed-node calculation, however, the energy is hardly improved by adding 
the Gutzwiller factor to the trial wave function. Considering also other calculations we 
performed, we feel that this is a rather general feature of the fixed-node method. In the 
fixed-node method, the energy that is obtained is determined by the sign structure of the 
trial wave function, via the ratios of the wave function in neighboring configurations with 
different signs. In order to significantly improve the energy, one must be able to change the 
sign of the trial wave function in individual configurations, as that is the most important 
factor determining the energy that can be obtained. By a Gutzwiller factor (or, similarly, 
a Jastrow factor) one does not change the sign, but only the ratios of the wave function. 
Another interesting point is that the Monte Carlo calculations become less stable for larger 
U. This is due to the fact that the parameter r has to be taken smaller in that case, such 
that convergence is slowed down. For the fixed-node method this is not a very big problem, 
as one can simply sample for a longer time to obtain a result with desired accuracy, but in 
the nodal-relaxation method one can not afford to have longer paths in the calculations, as 
the sign problem starts showing up after a certain number of steps. 

Another problem concerns the observables that can be calculated. For operators that 
commute with the Hamiltonian, similar formulas as have been given for the energy in Sec- 
tion m can be used, although (in fixed node) the results obtained for operators other than 
the Hamiltonian itself can not present a bound to the exact result. For operators that do not 
commute with the Hamiltonian, however, there is a more serious problem; in that case, the 
mixed estimate does not yield the correct value of the observable for the ground state of the 
(effective or true, depending on the method) Hamiltonian. Usually, one uses an extrapolated 
estimate in that case, taking two times the mixed estimate minus the variational result 
(i.e., the value of the observable in the trial state). In practice, this leads to rather good 
results; however, one is unable to check how good the approximation is. 

Recently, Zhang et a/.0 presented the Constrained Path Monte Carlo method for lattice 
fermions, which they claim to be a better method than fixed node. Like we, they apply 
some kind of fixed-node principle, but in a continuous space of Slater determinants, which 
seems to cause that the CPMC energy is a better approximation of the true energy than 
the fixed-node energy. Also, they claim that they have a better method {backtracing) to 
perform calculations for non-commuting operators. It is not very clear to us how CPMC 
compares to fixed node combined with nodal relaxation. It seems to us that our method is 
more generally applicable, as we can use any wave function as trial wave function (as long 
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as it can be calculated relatively easily) and we can also treat systems of frustrated bosons 
or spins, while their method restricts the trial wave function to be a Slater determinant. 
In any case, it would be very interesting to make a better comparison of both methods for 
some different systems. 

Finally, we state that the most important ingredient of our method is the trial wave 
function. The quality of the trial function, i.e., its sign structure, determines how good the 
approximation is that can be obtained. This is not a trivial problem, for, as Ceperley states, 
' . . . little progress [has been made] in stating exact conditions that nodes must obey. Most 
studies of wave functions aim at improving the energy of the wave function itself, while for 
our purposes it would be necessary to explore its sign structure, in order the improve the 
energy that results after a fixed-node Monte Carlo simulation. We still not have a clear 
understanding, yet, of the requirements that a wave function must fulfill to make a good 
fixed-node trial wave function. However, we feel that the fixed-node method combined with 
nodal relaxation offers great new opportunities to tackle problems of lattice Hamiltonians. 
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TABLES 



TABLE L Comparison of the exact ground-state energy to the energy of mean-field and 
Gutzwiller wave functions, and to the energy obtained by MC simulations, for the 4x4 square 
(PBC), with 5 spins up and 5 down, for [/ = 4. A mean-field (MF) homogeneous (H) wave 
function has been used. The Gutzwiller wave function (GMF) was used with g = 0.6; its energy 
was calculated by variational Monte Carlo. For the fixed-node (FN) and nodal relaxation (NR) 
Monte Carlo calculations, a trial wave function has been used as indicated. Energies are given per 
site and in units of t. Constrained Path Monte Carlo (CPMC) results were taken from Ref. 11, 
exact results from Ref. |l^. 



u 


MF 


GMF 


FN/MF 


FN/GMF 


NR/GMF 


CPMC 


exact 


4 


-1.109 


-1.212 


-1.2186(4) 


-1.2201(4) 


-1.2234(6) 


-1.2238(6) 


-1.2238 



TABLE II. Various mean-field and quantum Monte Carlo calculations (QMC) of the exact 
ground-state energy of a 10 x 10 square (PBC), at half filling and with zero total 5"^, for U = A. 
Slave Boson MF result has been taken from Ref. O, QMC results from Refs. 14 and f^. 



method 




energy per site 


Mean Field (AF) 


-0.797 


Slave Boson MF 


-0.839 


Gutzwiller VMC 


-0.842 


QMC (Hirsch '85) 


-0.88(3) 


QMC (White '89) 


-0.864(1) 


FNMC 


-0.852(2) 


Power MC 


-0.860(5) 
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FIGURES 



FIG. 1. First part of a FNMC simulation, for a 2 x 2 x 2 cube (see text for details). Each 
star indicates the average of the energy samples in 10 Monte Carlo steps. The drawn line indicates 
the exact ground-state energy of the fixed-node effective Hamiltonian. The square at the right 
indicates the resulting energy, obtained after averaging over several thousands of steps, with a 
statistical error which is smaller than the symbol size. 

FIG. 2. Nodal relaxation for the system as in Fig. Two calculations using different trial 
wave functions are shown (see text for details). Error bars indicate the estimated values of the 
ground-state energy after each step of the calculation. 
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